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Energy minimization of Ising spin-glasses has played a central role in statistical and solid-state physics, 
facilitating studies of phase transitions and magnetism. Recent proposals suggest using Ising spin-glasses for 
non-traditional computing as a way to harness the nature's ability to find min-energy configurations, and to take 
advantage of quantum tunneling to boost combinatorial optimization. Laboratory demonstrations have been 
unconvincing so far and lack a non-quantum baseline for definitive comparisons. In this work we (i) design 
and evaluate new computational techniques to simulate natural energy minimization in spin glasses and (ii) 
explore their application to study design alternatives in quantum adiabatic computers. Unlike previous work, 
our algorithms are not limited to planar Ising topologies. In one CPU-day, our branch-and-bound algorithm 
finds ground states on 100 spins, while our local search approximates ground states on 1, 000, 000 spins. We 
use this computational tool as a simulator to study the significance of hyper-couplings in the context of recently 
implemented adiabatic quantum computers. 



I. INTRODUCTION 

The Ising spin-glass model was first proposed by E. Ising 
in 1925 as a mathematical model to understand the dynam- 
ics of phase transitions in ferromagnetic systems. Such sys- 
tems are composed of particles that can be in either of two 
possible energy spin states. These spins interact in pairs to 
produce an energy landscape that describes the overall behav- 
ior of the system. The model is described in graph-theoretic 
terms by representing atoms in a crystal with vertices and 
bonds between atoms with edges. Despite its simplicity, the 
model has become an essential research tool in the analysis 
of different kinds of physical systems such as stiff polymers 
ifTHl and genome sequences 0| which can be mapped exactly 
or approximately to an Ising model. Since physical systems 
found in nature are often disordered, unless cooled to 0°K, 
the model incorporates randomness — either in the realization 
of the atomic couplings or the spin states of the atoms. 

Definition 1 Spin glasses are solid materials in which the 
magnetic moments (spins) have disordered orientations and 
the strength of the nearest-neighbor spin interactions (bonds) 
are randomly distributed. 

Physical and chemical properties of a crystal depend on the 
total energy of the bonds, which depend on atomic states. In 
particular, the lower the total energy, the harder the material. 
Estimating total energy by a graph-based function facilitates 
the use of graph algorithms to study properties of solids. For 
instance, if an Ising graph represents a ferromagnetic system, 
then finding the generating function of cuts in said graph is 
equivalent to calculating the distribution of physical states 
over all possible energy levels. On the other hand, finding 
the minimum-cut (max-flow) of an Ising-model graph that 
represents an amino-acid sequence is equivalent to calculating 
the lowest-energy configuration for a corresponding protein 
0]]. For most Ising models commonly studied in the Physics 
community, the latter problem is equivalent to finding the 
ground-state energy of the underlying physical system. 
Formally, the ground-state determination problem (GSD) 



is defined as follows. Given an instance of an Ising-model 
graph, find the set of spin-state values or spin configuration 
that minimizes the overall energy of the underlying physical 
system described by the graph. Such a state is known as the 
ground state of the system. Barahona pj] proved that, for a 
general random-field model, the GSD problem is NP-hard. 
Thus, all known algorithms for finding optimal solutions 
exhibit super-polynomial runtime. 

Computing based on energy minimization in physical sys- 
tems. Given that many physical systems have a natural ability 
to find least-energy states, researchers are currently attempt- 
ing to exploit this phenomenon to perform useful computa- 
tion. At the atomic scale, in addition to high bit-density, en- 
ergy optimization can be aided by quantum tunneling, which 
effectively reduces the number of local minima. Thus, GSD 
problems are of particular interest to quantum-information re- 
searchers because they are suitable candidates for evaluating 
the performance of adiabatic quantum computers (AQCs). 
Recently developed AQCs employ an architecture based on 
Ising spin systems lfl9ll . First, the spin system is configured 
to represent a given combinatorial problem, i.e., the spin in- 
teractions are carefully controlled rather than random as in 
spin glasses. The ground state is found via quantum anneal- 
ing (the quantum analogue of thermal annealing), then read 
off as a bit sequence and interpreted as an answer to the com- 
binatorial problem. While the classical formulation of GSD is 
NP-complete, Oliveira and Terhal lfl2ll proved that formulat- 
ing a general GSD instance in the context of AQC is QMA- 
complete j22ll . Since the complexity of GSD is universal with 
respect to both quantum and conventional forms of computa- 
tion, it is important to consider how well an approximation 
to the ground-state energy can be obtained by purely clas- 
sical combinatorial optimization techniques. Consequently, 
Bansal et al. J2l proposed an approximation algorithm for 
GSD on Ising spin lattices, which essentially simulates these 
AQC architectures U(\ 17^/ . and thus limits their potential for 
quantum speed-ups. To approximate the least energy with e 
accuracy, the algorithm from Q] requires runtime exponen- 
tial in 1/e, which is hardly practical. In contrast, we propose 
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a branch-and-bound algorithm and a high-performance local 
search that quickly finds near-optimal energy values for arbi- 
trary Ising topologies. Such techniques can be used to study 
properties of solid-state materials, as well as critically assess 
the performance of non-traditional computing devices based 
on energy minimization in Ising spin-glasses. The main con- 
tributions of our work are summarized as follows. 

• A branch-and-bound algorithm for solving GSD exactly 
on Ising systems with up to 100 spins. 

• A high-performance local search algorithm for Ising 
spin-glasses. Empirical results show that it scales better 
than other GSD algorithms and produces near-optimal 
solutions for small- to medium-sized instances. 

• A generalization of GSD for simulating energy mini- 
mization in physical systems. In particular, we propose 
a self-contained number-factoring algorithm based on 
this approach. These results can be used as a baseline 
for evaluating the performance of non-traditional com- 
puting devices that solve hard problems via energy min- 
imization (e.g. AQC). 

• A comparison of two potential spin-glass architectures 
for AQC number factoring. 

The rest of the paper is structured as follows. In Section ITT1 
we discuss common variants of Ising models, as well as the 
best known algorithms for solving GSD. Sections UTTl and HVl 
introduce our algorithms for finding ground states. Section 
[V] reports empirical results for calculating ground states. We 
build upon these findings and describe our generalization of 
Ising models using hypergraphs in Section [VI] and compare 
two potential architectures for AQC integer factoring. We fi- 
nalize the discussion with concluding remarks in Section |yTH 

H. BACKGROUND AND PREVIOUS WORK 

Due to its flexibility, the Ising model has been reinterpreted 
by researchers to analyze different kinds of physical systems. 
Potts Il4ll suggested a generalization of the model where the 
spin values are uniformly distributed about the unit circle. 
In the Edwards-Anderson (EA) ^ interpretation, the spins 
are binary ±1 values, and the strengths of the atomic cou- 
plings are independent, identically distributed random vari- 
ables. Ising systems that share this property are known as spin 
glasses since the random positive/negative edge-weights sim- 
ulate the solid-state structure of chemical glass. Let Gi S i ng = 
(V, E) denote a spin-glass graph with n vertices. Each ver- 
tex ii 6 y is annotated with spin value S u £ {±1} and is 
assigned a magnetization weight h u . For u, v £ V, define 
(u, v) £ E to be an edge representing a bond between two 
adjacent spins with assigned weight J uv chosen randomly 
from either the standard Gaussian (p = 0,8 = 1) or the ±1- 
bimodal distributions. The internal energy of the system for a 
particular configuration of spin values a = {Si} is given by 

n 

E(a) = - J2 JijSiSj ~J2 hiSi (1) 

(i,j)€E i 



where the summation considers all pairs of adjacent spins. 
Putting together the energies of all spin configurations gives 
the Hamiltonian of the system. Thus, the ground state is given 
by E gs = mm(E(cr) | V a £ 7r„), where 7r„ is the set of all 
possible ra-spin configurations. Whether we are interested in 
the lowest-energy value or the n-spin configuration with such 
energy, | 7r n | = 2™ because each of the spins can take on one 
of two possible values. As discussed in the next section, en- 
ergy minimization is typically NP-hard. Thus, calculating the 
ground state exactly using an exhaustive search algorithm is 
feasible only for small Ising models. To provide a scalable 
way of finding ground states or approximating their energies, 
we need to employ heuristics such as those proposed in Sec- 
tion [TV] First, consider the following definitions. 

Definition 2 A bond is satisfied if and only if the config- 
uration of its incident spins minimizes its weight (coupling 
strength) such that —JijSiSj = — | J,j |; otherwise, the 
bond is unsatisfied. 

Definition 3 A set of spins S C V is frustrated if there is 
no configuration of the spins that satisfies all the bonds (u, v) 
connecting the spins in the set, i.e., v,u G S. 

Alternatively, the ground state is defined by the spin configu- 
ration that minimizes frustration in the Ising system, provided 
that the system is not affected by an external magnetic field. 
That is, in a spin configuration without any frustrated spin- 
sets, all bonds have been satisfied and the energy of the sys- 
tem is E gs = — Y2u j)£E I JiJ I' w hich is a lower bound of 
the energy function. 

In general, Ising-model graphs are not limited to a par- 
ticular topology, but two- and three-dimensional lattices 
are most commonly considered in the literature. To simu- 
late the behavior of infinite spin glasses, it is common to 
require that the spins lying on the dimensional boundary 
be connected to the spins on the opposing boundary on the 
same dimension. This can be viewed as a type of (peri- 
odic) boundary condition. In particular, only one periodic 
boundary condition is imposed for each dimension of the 
lattice. However, it is sometimes desirable to have boundary 
conditions on some but not all of the lattice dimensions. For 
example, a 2-D lattice may have zero (planar grid), one or two 
boundary conditions. When no boundary conditions are im- 
posed, some (boundary) spins have fewer than four neighbors. 

Complexity of GSD. In his work, Barahona [3fl showed that 
the NP-complete task of finding a maximum set of indepen- 
dent edges (edges with no common incident vertices) in a 
graph can be reduced to GSD on a cubic grid. Although 
most variations of GSD are known to be NP-hard, there are 
a few cases where the structure of the graph can be exploited 
to solve the problem in polynomial time. For example, Bieche 
et al. |6Q proved that the GSD problem on planar graphs 
with zero magnetization can be solved in polynomial time by 
showing a reduction to the minimum- weight perfect matching 
(MWPM) problem. It follows from their work that GSD in- 
stances with zero magnetization (hi = 0) and 0- or 1-periodic 
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TABLE I: Ising spin-glass properties that make the GSD problem NP-hard on lattices. MWPM stands for 
minimum-weight perfect matching. 



boundary conditions can be solved in 0(n 3 ) time 0). Specif- 
ically, the algorithm identifies the frustrated faces in a grid 
(4-cycles that have an odd number of negative edges) as ver- 
tices in a new graph Gf — (F,Ep). Gf is complete and 
each edge e = fj) 6 Ep is assigned a weight equal to the 
sum of the absolute weights of the edges in the original graph 
that are crossed by the minimum path that connects fi and 
fj. Recall that the edges in Gf are connecting sets of frus- 
trated spins. Therefore, minimizing the sum of the weights 
connecting fi £ F implies that we are minimizing frustration 
(Definition 3). Thus, finding the ground state is reduced to 
finding a MWPM on Gf- However, although MWPM is solv- 
able in polynomial time, the runtime is impractical for large 
instances and suffers from a big memory footprint due to the 
size of Gf- To overcome these limitations, the work in lfl3ll 
describes a heuristic based on the MWPM reduction where a 
reduced graph Gf is used instead of a complete one at the 
cost of sub-optimality. 

Table U shows the Ising lattice properties that make a GSD 
problem poly-time solvable and identifies the algorithms that 
are commonly used. Note that the number of dimensions, 
the number of boundary conditions and the presence of an 
external magnetic field are the main factors in determining 
whether an instance is NP-hard or not. More precisely, when 
we consider lattices with more than two dimensions or with 
two boundary conditions, the graph is no longer planar and 
the reduction to MWPM breaks down. Another special case 
considered in the literature is that of ferromagnetic ( Jjj > 0) 
GSD instances. Barahona [4] reduced this particular problem 
to (s-i)-min-cut or max-flow. 



III. FINDING EXACT GROUND STATES 

To better control the trade-offs between runtime and solu- 
tion quality of heuristics, it is important to design algorithms 
that are guaranteed to find exact ground states on smaller 
instances. The solutions obtained from such instances are 
used to evaluate scalable heuristics. 

Branch-and-bound (B&B). For general optimization prob- 
lems, B&B considers incomplete or partial solutions, where 
only a subset of the problem variables are assigned admissible 
values. Partial solutions are systematically constructed via 



branching. The branching process only develops partial 
solutions that are deemed promising, i.e., those that may lead 
to the optimal solution. Conversely, partial solutions whose 
cost is too high, are "bounded away" or pruned. 

B&B on Ising systems. Our B&B algorithm proceeds as fol- 
lows. First, all spins are labeled as unassigned-their value can 
be set in the future to either 1 or —1. The algorithm then cal- 
culates the lower bound of Equation|2] 

n 

E lb = - i J i,j \-J2\ hi i (2) 

{i,j)eE i 

It then selects a spin i and branches on one of the possible 
values for the spin that has not been explored yet. In each 
branch, the incremental change in En, caused by the assign- 
ment is recorded as follows. For each spin j adjacent to i that 
has already been assigned, increase (decrease) En, by twice 
the amount of the positive (negative) bond connecting i and j 
if they have opposing (aligned) spin values, 

!2 J i,3 S i S j if S i + S 3 and J i,3 > °' 

-2 J2 AjSiSj if Si = S and < 

Furthermore, the corresponding change due to the magnetiza- 
tion of the spin is also recorded, 

E s = [ if Si = -1 and h, > 0, 

lb \-2hi if Si = +1 and hi < 

Once the spin is assigned, the algorithm branches out to an- 
other spin and performs the same procedure. When all spins 
have been assigned, En, represents the energy of the spin 
configuration generated by the branching process. To con- 
tinue searching the configuration space, the branching pro- 
cess backtracks to the last assigned spin, flips the spin's value 
and updates En,. If both spin values have already been tried, 
then the algorithm continues backtracking while relabeling 
the spins as unassigned. Since each spin can take one of two 
values, this branching process generates a full binary search 
tree where the leaves correspond to all possible spin configu- 
rations in the Ising system. 
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Initially, we use a linear-time greedy approximation (E gs ) 
as our bounding value. During the branching process, if the 
energy of the partial solution exceeds E gs , then we can safely 
prune this branch and backtrack to the previous assigned spin 
without making any further assignments. The algorithm either 
tries the opposite spin value or backtracks again if both spin 
values have already been tried. If the search assigns all the 
spins in the graph and the corresponding minimal energy state 
is lower than E gs , then we set E gs to this new energy value. 
After searching all promising branches, E gs will assume the 
ground-state energy. This standard bounding technique alone 
improves the scalability of the branching process by an order 
of magnitude over exhaustive search (see FigureQ]). 

To further improve the scalability of our B&B algorithm, 
we designed a prune-by-dominance technique that consists 
of identifying partial solutions whose partial energy can be 
improved (lowered) by modifying the configuration of cur- 
rently assigned spins. Note that, whenever we assign a spin s, 
there is a set F s of spins adjacent to s for which all neighbor- 
ing spins (including s) have also been assigned. Early in the 
branching process, F s is likely to be empty since only a few 
spins have been assigned. As spins are assigned, the set F s 
increases. Figure [2] shows an example of s and F s on a small 
grid. The size of F s is no greater than the degree of s. Note 
that the energy of the spins in F s is localized in the sense that 
it will not be affected by further spin assignments. 

Lemma III.l Let F s be the set of spins such that Vi G F s , all 
spins adjacent to i are assigned. Then the partial energy of the 
spins in F s will not be affected by additional spin assignments. 

Proof: Let i G F s , then the partial energy lower bound 
that is localized around spin i is the given by Ej b = 
^J2(i j)£E JijSiSj ^ 2/ij (the ± stands for the cases de- 
scribed in Equations [3] and 0). By definition of F s , we know 
that every spin j adjacent to i has also been assigned. Sup- 
pose that later in the branching process we assign spin u and 
this causes a change in E\ h . By definition of E\ h , u must be 
adjacent to i. This implies that i ^ F s since its neighbor u had 
not been assigned until recently, which is a contradiction. 
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This allows us to evaluate a partial solution by flipping the 
values of the spins in F s and comparing the partial energies. If 
any of the partial energies are lower than En, or if the modified 
configuration corresponds to a visited partial solution, then we 
know that the current partial solution is unpromising and we 
can safely backtrack. 



FIG. 2: Illustration of s and _F S on a small grid. 

Lemma III.2 Let F s be the set of spins such thatVi G F s , all 
spins adjacent to i are assigned. Also, let <7f 3 correspond to 
some configuration of the spins in F s . If we can find a' F such 
that E(a' F ) < E(o~f 3 ), then any branches that include o~f 3 
can be safely pruned. 

Proof: Consider the partial energy E(aF 3 ) and the total 
energy of any complete spin-configuration a that extends the 
partial spin-configuration. Let a' F be a configuration of F s 
that minimizes the partial energy, i. e., E(a' F ) < E(<jf 3 ). By 
Lemma llII.il new spin assignments will not affect E(<jf b ). 
Since <jp s is included in a we can swap ap 3 with a F so 
that E(a) is also minimized. Thus, any partial or complete 
spin-configuration that includes <tf s is not promising. 

Observe that in cases when different branches are unlikely 
to have equal partial cost (e.g., when couplings and magne- 
tizations are random), for two branches, the probability that 
the first branch dominates the second branch is approximately 
1/2. Let < c < 1 be the fraction of (2 k partial solu- 
tions) that require branching. Then we can expect to prune 
c(2 fe /2) = c(2 k - 1 ) of these branches. As seen in Figure Q] 
this pruning technique improved the scalability of our B&B 
algorithm by 1-2 orders of magnitude, allowing it to solve 
100-spin lattices in a day. However, even with the techniques 
proposed, B&B takes exponential time in the worst case and 
therefore fails to scale beyond 100 spins. 
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FIG. 1: Performance of B&B techniques on 2-D spin lattices. 



IV. GSD THROUGH LOCAL SEARCH 

Due to the difficulty of solving general instances of GSD, 
many researchers |@] lfl3ll have developed heuristic methods 
to improve the scalability of their algorithms at the expense of 
solution quality, typically based on slow Monte Carlo simula- 
tions. However, because of the role that Ising models play in 
simulating real-world phenomena, it is desirable to have much 
faster techniques. To this end, we propose a high-performance 
local search that meets such scalability and performance re- 
quirements. 

Our local search is an iterative improvement algorithm 
that modifies the bipartition induced by an arbitrary spin 
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configuration (positive spins are placed in one partition 
and negative spins in the other). The algorithm performs a 
sequence of incremental changes to the bipartition, organized 
as passes. These changes consist of spin moves that place 
a particular spin in the partition opposite to the one it is 
currently in. At the beginning of each pass, the energy 
differential (gain) of performing each possible move is 
calculated. A positive gain implies that the move decreases 
the overall energy while a negative gain increases it. During a 
pass, the move that produces the largest gain is selected and 
executed. The corresponding spin is then labeled as locked, 
i.e., it cannot be selected again in the current pass to prevent 
"undo" moves. The pass continues selecting and executing 
the best moves until all spins have been locked. At the end 
of the pass, we save the best-seen bipartition produced by 
the sequence of moves. This bipartition is then used as 
the starting solution of the next pass. The entire algorithm 
terminates when a pass fails to obtain an improvement in 
energy as shown in Figure [3] Note that, in the absence of 
positive-gain moves, a negative-gain move can be selected. 
Thus, a pass may accept a solution that is worse than the 
existing solution (hill-climbing). This helps to reduce the 
probability of getting trapped in local minima. Figure [5] 
illustrates the progress of our local search in terms of solution 
costs during individual passes on a 1024-spin glass. The 
initial random solution used in our local search is generated 
in linear time. Assuming that initial spin configurations are 
drawn from a uniform distribution, our local search finds an 
optimal solution with probability at least 1/2" for n spins (in 
practice, much higher than that, as indicated in Figure|9]l. The 
speed of our algorithm can be converted into better solutions 
by generating independent random inital spin configurations, 
running (otherwise deterministic) optimization passes on 
each, and selecting the best result. 

Efficient gain updates. Each move causes a change in the 
local energy surrounding the selected spin, therefore, the 
gains of the neighboring spins need to updated after each 
move. When the graph capturing a spin system is sparse 
(e.g., Ising lattices), only a constant number of gain updates 
are executed per move. These gain updates are performed 



Input: Ising spin-glass graph G iaing 
Output: Approximate ground-state energy E* 



Spin Partition SP* := RAND.SPINJ>ARTITION{G iaing ) 
while solution quality improves do 

Gains container GC ■= COMPUTE.GAINS(SP*) 
Spin Partition SP := SP* 
while GC has unlocked spins do 

Move m := SELECT J3EST.MOV E(GC) 
APPLY. MOV E(SP, m) 
UPDATE.GAINS(GC, m) 
LOCK.SPIN(GC, m) 
if energy decreased then SP* := SP 
end while 
end while 

return ENERGY(G iamg , SP*) 



FIG. 3: Pass-based local search with hill-climbing 



FIG. 4: Heap-based data structure for performing gain updates ef- 
ficiently. Dashed lines indicate pointers to child nodes. Spin array 
pointers are updated accordingly after every heap swap. 

efficiently using the custom heap-based data structure shown 
in Figure |4] The data structure consists of two arrays. The 
first array implements a traditional binary heap while the 
second array allows quick access to the heap-array element 
that contains the gain-update value of a particular spin. To 
perform gain updates, we can access the specific value in 
O(l) time, update the value, and perform the necessary swaps 
to maintain heap-order property. Since only log(n) swaps 
are required in the worst case (where n is the number of 
spins), our data structure allows us to perform gain updates 
much faster than naive implementations that require scanning 
the entire set of n gain values. Since a total of n moves 
are performed during a pass. This gives a total runtime of 
0(n log(n)) for a single pass. 



V. EMPIRICAL VALIDATION 

We evaluated single-threaded implementations of our 
algorithms on a conventional Linux server, although our 
local search is trivial to parallelize to a multicore system or a 
distributed cluster. For 15x15 spin lattices our local search 
finds exact ground states in 95% of independent random starts 




-2500 u 1 1 1 1 ■— 1 

200 400 600 800 1000 

Moves performed in a given pass 

FIG. 5: Progress of local search in terms of total energy during indi- 
vidual passes on a 1024-spin glass. The lowest-energy state observed 
during each pass (Ex = -2456, E 2 = -2463, E 3 = -2465, re- 
spectively) becomes the starting state of the next pass. 
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(exact solutions were obtained from |21[|), otherwise solutions 
are 5% sub-optimal on average. Figure |7Ja) compares the 
average solution quality of local search for 2-D spin glasses 
with Gaussian-distributed couplings and hi = (instances 
with hi ^ are not allowed in 12110 . For each instance we 
considered four different levels of effort with an increasing 
number of independent random starts. To obtain the average 
solution quality we computed 1000 output samples using 
1, In 2 n and n random starts (n is the number of spins) per 
instance. For n 2 random starts, we used fewer output samples 
and provide confidence intervals. As expected, the solution 
quality improves as the number of random starts increases. 
When at least In 2 n random starts are used, our heuristic 
produces high-quality solutions (> 95%) for five of the 
benchmarks while its runtime does not exceed 17 seconds for 
the largest benchmark (2500 spins). Note that the expected 
solution quality slowly decreases for larger instances. Figure 
|7|b) shows similar results for benchmarks with ±l-bimodal 
coupling distributions, but solutions are closer to optimal. For 
all but one of the benchmarks, our heuristic requires only a 
single random start to find high-quality solutions. 

Local optimality. We verified that the configurations re- 
turned by our heuristic cannot be improved by modifying a 
small number of spins. We used our B&B to find optimal 
configurations of groups of 25-49 adjacent spins within larger 
configurations. In our experiments, the solutions produced by 
our local search were never improved by this technique. 

Runtime. Figure [6] shows that our heuristic scales to a mil- 
lion spins and its runtime is consistent with the complex- 
ity estimate of 0(n\og(n)) per pass (see Section ITVli. We 
compared the runtime of our local search against that of 
the MWPM-based heuristic proposed in lfl3ll (see Section 
Hil l. Recall that the heuristic works on a reduced dual graph 
instead of the complete one required to find exact ground 
states. The reduced dual graph ignores those edges that are 
deemed "heavy", i.e., with weights above a chosen thresh- 
old Cmax — c ■ Jmax, where c = 2,4,6. The idea is that 
these heavy edges are rarely contained in the optimum so- 
lution and can be ignored. Thus, the solution quality of the 
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FIG. 6: Runtime of local search on large 2-D spin lattices. 



heuristic depends on c max and the positive-to-negative edge 
ratio since the more edges we ignore, the higher the proba- 
bility of missing the ground state. In contrast, our heuristic 
does not have such a dependency and works on all instances. 
Furthermore, it is not clear how to chose c max in the case that 
the couplings weights are Gaussian distributed. In contrast, 
our heuristic can be used with any coupling distribution. Ta- 
ble IV in lfl3ll shows the runtime and solution quality of the 
MWPM-based heuristic for different values of c max on ±1 
grid graphs of size 164 x 164 with .5 positive-to-negative edge 
ratio. The fastest runtime of the cited heuristic is obtained 
when c max = 3 taking an average of 58.72 seconds (with 
negligible deviation) and producing the optimal value only 
61% of the time. By comparison, our local search heuristic 
on a comparable benchmark with a single random start takes 
about 8.5 seconds. Thus, we can perform 7 random starts in 
the same period of time. However, we have no way of com- 
paring the quality of our solutions since we do not have ac- 
cess to the same benchmarks. Section 4.5 of ]9j] describes 
the branch-and-cut approach used by the Cologne Spin Glass 
Server. This algorithm is limited to lattices without magneti- 
zation (our techniques do not suffer this limitation) and men- 
tions a 128-second runtime for the largest instance (50 x 50). 
In contrast, our local search heuristic takes only an average 
of 16.63 seconds using ln 2 (n) random starts on the same in- 
stance and produces solutions that are within 95.5% of the 
optimal value on average. Solutions closer to optimum may 
be found if the number of random starts is increased. 



VI. GSD FOR ARBITRARY HYPERGRAPHS 

As discussed in Section |U recent AQC architectures lfl5ll 
are based on Ising spin glasses. Implementations described by 
D-Wave S yste ms use non-planar topologies, and recent exper- 
iments in 111811 demonstrate direct coupling of more than two 
spins. Hence, we extend the conventional spin-glass model to 
use hyper-couplings that connect a set of at least two spins. 
The new energy function is given by 

E ( a ) = - 2 je n s * - J2 kiSi (5) 

This formulation further extends the applicability of our 
computational approach to energy minimization. Adapting 
our algorithms to handle this formulation allows us to study 
a greater variety of physical systems and solve a wider range 
of combinatorial problems (e.g., number factoring-discussed 
later in this section). Our algorithms required minor modifi- 
cations to handle hyper-Ising models. In the case of B&B, 
the incremental energy calculations (Equations [3] and |4|i 
were updated to conform with Equation [6] For local search, 
the gain-update procedure was revised to account for the 
presence of hyper-couplings. 

Number-factoring as an optimization problem. We observe 
that integer factorization is equivalent to the optimization of 
the following function, 

f(x,y) = (N-xyf (6) 
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FIG. 7: Expected solution quality (100% corresponds to the exact solution obtained from |21jP 
produced by using different numbers of random starts on 2-D spin lattices. 



where N is the odd integer we wish to factor and x, y are pos- 
itive integers. The minimum of f(x, y) is reached when x and 
y are factors of N. To solve this optimization problem using 
GSD, we must construct an Ising system whose ground state 
corresponds to the global minimum of f(x, y). Since spins 
are ±1 binary variables, we re-formulate f(x, y) in terms of 
binary digits. Let n x and n y be the number of binary digits 
required to represent x and y, respectively. 



f(x,y) 



N 




(7) 



Since N is odd, x and y are also odd and therefore the last 
term in each sum is one. Thus, the variables xq and yo are 
ignored. The number of spins is (n x — 1) + (n y — 1) =| V \. 
Setting S k = 2x k - 1 and S, lx+k = 2y k - 1, 



f(x,y) 



N 



2 1 Sn ~ n * +i 



2 « y 



-1 1 $n—n x 



1-Si 

2 r-^ + 1 



(8) 



The magnetization weights and coupling strengths are given 
by the product expansion of Equation [8] Terms of the form 
c • Si correspond to spins with magnetization hi = c, and 
terms of the form c • SiSj..S k correspond to hyper-couplings 
e = (i, fc) with strength J e = c. Figure [8] shows the Ising 
graph constructed for factoring 21. 



number N, which is also typical behavior in adiabatic 
quantum computers (AQC). Figure [9] shows the output 
probability distribution of factoring the number 612171, 
which has several factors and is therefore easy to factor. The 
probability of success drops sharply for semi-primes, e.g., 
the probability of factoring 580003 in one attempt is only 
.005. Therefore, our heuristic takes longer to find correct 
factors in such cases. The current implementation of our 
algorithm is not particularly optimized, but factored the 
semi-prime 10185081163 = 100511 x 101333 using 15,000 
random starts in 13 seconds. While further optimizations 
may significantly improve runtime, present results provide 
a computational baseline for performance evaluation of 
non-traditional computing devices that solve hard problems 
using the principle of energy minimization. Such devices, 
e.g., AQC, would need to improve on our results by producing 
distributions that increase the probability of factoring the 
correct number. 

Significance of hyper-couplings. Recently, Peng et al. lfl5ll 
implemented an AQC number-factoring algorithm in NMR 
technology using a Hamiltonian that ignores hyper-couplings. 
By simulating the functionality of the same algorithm, we 



Computational experiments. B&B can factor up to 21-bit 
numbers in approximately 5 minutes. While leading-edge 
number-factoring algorithms can do better, the results con- 
firm that B&B is general enough to work on hyper-Ising 
models. We tested our GSD heuristic by factoring specific 
numbers. In some runs, this technique produces the factors 
of N ± 2 or other incorrect numbers. Therefore, many 
independent random runs may be required to factor a given 



FIG. 8: Ising-model hypergraph for factoring N = 21. In this exam- 
ple, we assume x < y and set n x — 1 and n y = 2. The optimization 
function is f(x,y) = 210 + 84S + 88Si + US 2 - 20S Si - 
IOS0S2 + 2OS1S2 - WSoSiSi. 
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613407 



CO 
O 



0.08 
0.07 
0.06 
0.05 
0.04 
0.03 
0.02 
0.01 





578867 




OBUUU3 

Number factored 

(b)580003 is factored with probability .005. 
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FIG. 9: Probability distribution of factoring numbers using our GSD heuristic. The black lines show the probability of 
finding the correct factors using a single random start. Given that 100 attempts take < .01 seconds, such numbers can 
be factored reliably. 



replicated their experiment and explored factorization of 
larger numbers. As in lfl5ll . the number 21 was successfully 
factored even in the absence of hyper-couplings. However, 
our experiments show that, in general, omitting hyper- 
couplings alters the ground states so that correct factors 
cannot be found, e.g., for 35 and 91. Instead, some other 
numbers are factored such as 33 and 95. Although hyper- 
couplings have smaller individual weights than two-spin 
couplings, the number of hyper-couplings scales as 0(n 4 ), 
and their total weight eventually dominates f(x, y) for larger 
N, Control of three-spin hyper-couplings has recently been 
demonstrated J3, but only for adjacent particles, which 
would be insufficient for number-factoring applications. Fur- 
thermore, Equation |8] still requires four-spin hyper-couplings 
which, as current research suggests, are difficult to realize 
experimentally. 



olating the factorization equations. This penalty is minimized 
subject to the fixed values of m,\...m n . Equation[9]maps to a 
spin system with only two-spin couplings. Unlike the (uncon- 
strained) formulation with hyper-couplings, this formulation 
includes spins with fixed values and thus requires additional 
technology support (e.g., optical pumping of trapped ions). 
We used our algorithms to compare the use of hyper-couplings 
(assuming technological feasability) to the expansion from 
II 1611 . Figure[l0]compares output probability distributions gen- 
erated when factoring 51. The technique from lfl6ll produces 
a flatter distribution with a wider range. This is because the 
solution space is more complex and includes configurations 
with inconsistent ancilla spins (e.g., carry spin Cij = 0, but 



the partial product spin pi , 



-1,3 



* 1 = 1). Thus, this tech- 



Avoiding hyper-couplings via ancilla spins. The work in 
lfl6ll shows a method for expanding Equation [8] that avoids 
the use of hyper-couplings at the cost of a quadratic increase 
in spins. The new Hamiltonian is based on the factorization 
equations that are generated by decomposing long-hand bi- 
nary multiplication [lfj. Equation [8] is modified by intro- 
ducing ancilla binary variables. Let N = mim,2---m n . Let 
x = x\X2---Xk and y = yiy 2 ...y n -k- Then let p ?J denote 
the sum of pairwise products between the binary variables in 
Equation [8] and let Ci.j denote carry variables 



h n—k 
i=l j=l 

-Pi+i,j-i - 2q-i.j - 1/4) 2 - 1/8 



(9) 



The bits of N are linked to the bits of x and y by a series of 
equalities such as pi.o = mi and pk+i,j-i = rrik+j (see lfl6ll 
for details). f(x, y) now computes a penalty function for vi- 
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FIG. 10: Output probability distributions for factoring 
51 using the techniques from fT^l and fTr3l as simulated 
by our algorithms. 
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nique sometimes returns trivial decompositions of prime num- 
bers (e.g., 31), whereas direct use of hyper-couplings always 
results in proper factorizations. In summary, our experiments 
suggest that the expansion from lfl6ll is not computationally 
competitive with direct use of hyper-couplings. 

VII. CONCLUSIONS 

The problem of finding the ground state of spin glasses 
(GSD) has played a central role in statistical physics. Ad- 
ditionally, high-performance GSD algorithms can help evalu- 
ate architectural alternatives in adiabatic quantum computing. 
Recent work in 11711 concludes that quantum annealing loses 
to classical simulated annealing in a head-to-head compari- 
son, but can be improved by making different architectural 
choices. Although our algorithms do not simulate quantum 
optimization directly, they allow one to study problem reduc- 



tions and identify potential obstacles to successful optimiza- 
tion. The proposed B&B algorithm can find ground states for 
general 2-D spin lattices with up to 100 spins in 24 hours. For 
GSD instances with 1, 000, 000 spins, our local search heuris- 
tic obtains approximate solutions in < 4 hours. It provides a 
scalability advantage over conventional Monte Carlo methods 
and is not limited to special classes of GSD instances. When 
our heuristic does not find a ground state, it usually approxi- 
mates the least energy within 5%. 

We replicated a recently published empirical result in AQC- 
based number factoring 1151 . where the omission of spin-spin 
hyper-couplings did not undermine overall results. However, 
we have shown that, in general, omitting hyper-couplings will 
produce incorrect results. Furthermore, we demonstrated that 
techniques avoiding the use of hyper-couplings lfl6ll blur the 
output probability distribution, hamper finding correct factors, 
and require more repeated attempts to achieve success. 
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